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The dynamics of a ring of vortices in two-dimensional Bose-Einstein condensates (with and without 
an additional vortex at the center) is studied for (1) a uniform condensate in a rigid cylinder and (2) 
£N| a nonuniform trapped condensate in the Thomas- Fermi limit. The sequence of ground states (within 

these single-ring configurations) is determined as a function of the external rotation frequency by 
comparing the free energy of the various states. For each ground state, the Tkachenko-like excitations 
and the associated dynamical stability are analyzed. 

PACS numbers: 03.75.Kk, 67.40.Vs, 31.15.Pf 

r-| ■ I. INTRODUCTION 

The dynamics of quantized vortices has attracted great interest. Until recently most research has been done in the 
context of superfluid 4 He pj. The physics of an unbounded superfluid with an infinite number of vortices is relatively 
simple: In equilibrium the vortices form a two-dimensional rotating triangular vortex lattice, which supports small- 
amplitude collective modes 2j. In a bounded superfluid, however, determining the equilibrium configuration is not 
easy even for a small number of vortices. 

Havelock studied the dynamical stability of a ring with a small number of vortices. Subsequently, Hess 
studied the energy of various configurations analytically to determine the sequence of low-lying equilibrium states, 
and Campbell and Ziff p| evaluated the energies numerically for larger numbers of vortices. After heroic efforts, 
Yarmchuk et al. Q photographed some stable vortex patterns in superfluid 4 He. 

For a large number of vortices, Fetter developed a theory of the dynamics of a lattice of slightly displaced rectilin- 
ear vortices including three-dimensional distortions. Later, the description was extended to provide a continuum 
picture of the vortex dynamics |8( • Baym and Chandler formulated similar macroscopic hydrodynamics of rotating su- 
perfiuids . These various theoretical descriptions included the vortex-interaction effects discussed by Tkachenko @] • 
■ For a comprehensive review of vortex dynamics in rotating superfluids, see, for example, Ref. [To|. 

The strong interaction between helium atoms makes it difficult to manipulate the condensate and to observe vortices. 
Thus, the detection of quantized vortices in trapped Bose-Einstein condensates offers an opportunity to 

study the dynamics of superfluid vortices. In addition, the interaction parameter is tunable by varying the number 
of atoms, which helps clarify their many-body character and the detailed connection to fluid dynamics. Recent 
experiments at JILA 1141 measured the transverse elastic modes of the bounded vortex lattice that have been studied 
theoretically by d Hi El 13 . These Tkachenko modes in a large but finite lattice require appropriate and 
i-^J ■ nontrivial boundary conditions. 

Here, we study the equilibrium states and small oscillations of (1) a ring of N v vortices and (2) a ring of N v — 1 
O ■ vortices plus a central vortex for two specific cases: a cylindrical container with uniform fluid and a two-dimensional 
condensate in a harmonic radial trap in the Thomas- Fermi (TF) limit. These vortex states can be stable up to 
, N v ~ 9, beyond which other more complicated configurations become favorable. The associated normal modes are 
small-system analogs of the Tkachenko modes for a large vortex array. 
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II. TIME-DEPENDENT VARIATIONAL FORMALISM 



A macroscopic order parameter (the condensate wave function) = with associated particle density n — \^\ 2 , 
provides a simple description of a Bose-Einstein condensate near T = K. In the presence of an external trap potential 
Vt r , the condensate experiences both Vtr and the self-consistent Hartree interaction Va — g\^>\ 2 = gn arising from the 
interaction with all other particles, where g = 47rh 2 a s /M is the interparticle interaction parameter (assumed to be 
repulsive) and a s is the positive s-wave scattering length 0, |5(J • 

At zero temperature in the laboratory (nonrotating) frame, the macroscopic order parameter ^ satisfies the time- 
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rm= idv^U*^-^*), (4) 



dependent, nonlinear Schrodinger equation (Gross-Pitaevskii equation), 

^ = -|^ + K r * + sl*| 2 *. (1) 

If the external potential rotates at an angular velocity S~2, it is necessary to transform to a co- rotating frame, in which 
case the time-dependent Schrodinger equation becomes [2l| 

ih— = - — y + V tI *-n-rx P y + g\y\ 2 y, (2) 

where rxp = — ihr x V is the angular-momentum operator. 

For many purposes, however, it is more convenient to use an equivalent Lagrangian formalism based on the La- 
grangian functional |22l ]. 

£[#] =T [#]-£'[#]. (3) 

Here 

TK&l = / dV - , . 

2 V dt dt 

is the time-dependent part of the Lagrangian, the analog of the kinetic energy in classical mechanics. In a nonrotating 
frame, £[&] is the energy functional. In the reference frame that rotates at the same angular velocity O as the external 
trap potential, the appropriate free energy functional acquires an additional term, £'[$>] — £ — fi • L[\If], where 
L[^] is the angular momentum of the condensate. When l~i = flz is oriented along z, £' can be written as 

£'[*]= J dV ^\V^ 2 + V^\ 2 + 9 ^\ A + . (5) 

The action S associated with the above Lagrangian is the time integral S[^>] = J dt It is stationary with 

respect to small variations of iff and iff*, when iff satisfies the time-dependent Gross-Pitaevskii (GP) equation J5J. 
Thus the time-dependent GP equation is the Euler-Lagrange equation for the Lagrangian functional. 

This Lagrangian formalism has the following significant advantage: not only is this approach exact, but also it 
provides the basis for a powerful approximate variational method. If the condensate wave function depends on 
various parameters {p} of the given system, the resulting Lagrangian £({p}) provides the dynamical equations of 
motion of these parameters {p}. Moreover, when the system under consideration is in equilibrium and the variational 
parameters of the problem are stationary, the next-order perturbation (variation) of the Lagrangian yields the normal 
modes of the system in terms of these parameters. 

Here we use this method to study the dynamics of vortices in a ring in a two-dimensional Bose condensate. The 
positions of the vortices {ij} serve as the time-dependent variational parameters for two distinct cases: (1) a uniform 
condensate in a rigid cylinder and (2) a nonuniform trapped condensate in the TF limit |23|. 

In our variational approach, we use a trial wave function of the following form 

*(r, {r.,}) = |*(r)|e l5(r ' { ^ }) , (6) 

where \iff\ 2 is the density profile of either the uniform condensate or the nonuniform TF condensate. It satisfies the 
normalization condition 

N = J dV\V\ 2 , (7) 

where N is the number of condensate particles. In addition, ^(r, {r^}) is the phase of the condensate wave function, 
and the particle velocity field is given by v(r) = (Ti/M) VS'(r). The trial wave function © assumes that the vortices 
do not affect the density profile of the condensate (this approximation holds for small cores, which is valid in the TF 
limit); instead, the vortices directly determine the velocity field. For simplicity, we consider only a condensate in two 
dimensions, when the various terms in the Lagrangian are interpreted per unit length. 



III. FREE ENERGY AND EQUILIBRIUM STATES 

To proceed, it is first necessary to choose a specific form for the phase S of the condensate wave function. The 
boundary conditions are quite different for a uniform condensate and for a TF condensate, and it is convenient to 
treat them separately. 
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A. Uniform condensate in a cylinder 

The uniform condensate in a rigid circular cylinder of radius R must satisfy the boundary condition that the normal 
component of the particle current vanishes at the surface of the cylinder. This condition can be ensured with an image 
vortex of opposite circulation at the external position r'- = (i? 2 /r 2 ) Tj for each vortex at rj in the cylinder. 

The time-dependent part T of the Lagrangian can be easily evaluated. The trial condensate wave function is 

*(r,{r,}) = v^e <S(p,{r ' }) , (8) 

and 

%{'i})=E s i( r ' r i)+E 5 i ( r > > w 

where 

Sj(r, rj) = arctan ^] , (10) 

\ x — x j J 

S'(r,r>) = (-1) arctan (j^j , (") 

n = N/V is the particle density and N v is the number of the vortices in the condensate. In cylindrical coordinates, 
we obtain 

N v N v 

ni*j}) = E 7 ^) = E C 1 - ( 12 ) 

where each quantity is per particle and written in dimensionless units: the length scaled by the radius of the cylinder 
R, the time scaled by fl^ 1 = MR 2 /h and the energy by Tiff- 
in the free energy part of the Lagrangian, the quantity of interest is the difference AS' of the free energy between 
the vortex state and the vortex-free state. In our approximation, this difference arises from the supcrfluid velocity 
field v caused by the presence of the vortices. It can be simply determined as 

A£' = A£-QL Z , (13) 

where AS is the contribution from the superfluid velocity field in the nonrotating reference frame 

A£ = J d 2 r|*| 2 iv 2 , (14) 

and L z is the z component of the angular momentum 

L z = [ d 2 r |*| 2 z r x v. (15) 



In evaluating the free energy part, it is convenient to use a stream function \ instead of the phase of the condensate 
wave function S (which is effectively a velocity potential). Thus, we choose 

Ny Ny 

X(r,{r j }) = ^Zx j (r,r j ) + Y^ X ' j (r,r' j ), (16) 

where 

X,(r,r j ) = ln|r-r j |, (17) 
X ;-(r,r;-) = (-l)ln|r-r;.|, (18) 

and the velocity field is given by 



v(r) = z x Vx- 



(19) 
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This two-dimensional stream function has the advantage being single-valued. When describing irrotational flow in 
the presence of vortices, x satisfies Poisson's equation 



N v 



V 2 X (r,{r j })=2^6^ (r - r 3 -) . 



(20) 



The right-hand side vanishes everywhere except for the finite number of singular points at the location of the vortices 
(the image vortices lie outside the physical region and do not appear here). Thus, this stream function is essentially 
a superposition of two-dimensional Green's functions. 
Each individual term Xj ( r > r j ) has the Fourier expansion, 



Q(r — r,j 
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where 



(21) 



e(r-rj) 



if r > rj, 
if r < rj, 



(22) 



is the unit positive step function and 4> is the angle between r and Yj. Straightforward analysis then yields the free 
energy difference 



A£'({ rj }) =N v ln 



N„ N v 

£m(i-r 2 )-n£(i- 

3=1 J=l 



JV,, 



i<j 



1 — 2nrj cos 4>ij + rfr 



*7 



2J2 



(23) 



2r,r,- cos ( 



where ^ is the healing length of the condensate (the vortex-core-radius cutoff) and <pij = <f>i — 4>j is the angle between 
the ith and the jth vortices. Each vortex interacts with all other vortices including images, as is evident from the last 
term. Note that £1 is also written in dimensionless units scaled with Qq. 



B. TF condensate 



In the TF limit, it is more convenient to work with the equivalent grand canonical ensemble, introducing the 
chemical potential /x that is determined by the normalization condition N — J dV | \I/ 1 2 . Then the condensate wave 
function gives the following TF density profile 

|*(r)| 2 =n (l-r 2 ), (24) 

where uq = fj,/g is the central density and r is a dimensionless radius, scaled by the TF radius R = \J1[ijMuj\ (in 
the TF limit, \i ^> hcu±, which implies that R is large compared to the oscillator length y/h/Muj±). 

In the TF case, the condensate density vanishes at the boundary (r = 1). Thus the particle current automatically 
satisfies the boundary condition, and image vortices are unnecessary. When the energy is scaled by the characteristic 
energy HSIq — h 2 /MR 2 = (hu)j_/2fi) hujj_, the time-dependent part of the Lagrangian per particle is 

^{ri})=E^(ri) = -£k 2 
The free energy difference per particle can be evaluated using the same Fourier expansion (|21|l for the stream 
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functions (see Appendix lA"|l . which yields 



i N - r 

A£'({r,-}) = -£ (1 



i=i 



21n(|)+ln(l-r?) 



(2r?-l) -Ofl-rf) 2 



t<3 



E 



1 - 2nrj cos 0y + rf r| 



(26) 



r^r,- sin cpu arctan 

v 1 — r^rj cos < 



2 2> 

r — r 



l<0 



C. Euler-Lagrange equations of motion 



The Euler-Lagrange equations for the Lagrangian C = T — £' are 

d dC _ dC_ 

dt dfj drj ' 
d dC d£ 
dt BL dd>4 ' 



(27) 
(28) 



Since C is independent of tj for the configurations studied here, the first equation reduces to dC/drj — 0, which 
determines the precession frequency of the jth vortex [see Eqs. I|12() and l|25|l ]. 



fdAS 




(29) 



If the system is nondissipative and all vortices lie symmetrically spaced on a single ring of radius r, they are locked 
together and precess at the same rate <p = fair, ^)- 



D. Uniform precession of a single vortex 

It is interesting to start with the simplest case of a single vortex, which clarifies the distinct role of the factor In (R/£) 
between the Lagrangians for the uniform condensate and the TF condensate. Consider the uniform condensate, where 
the free energy part of the Lagrangian 112MI) contains the logarithm In (i?/£) only as an additive constant N v In (R/£), 
independent of the coordinates of the vortices. Thus, this logarithmic factor ln(i?/£) affects the energy-related 
properties, such as the ground-state configuration and the angular velocity for transitions between adjacent states; 
in contrast, the same factor does not affect either the radius of the ring of vortices or the precession frequency. For 
example, Eq. I|29|l gives the scaled precession frequency for a single vortex at a distance r\ in a uniform nonrotating 
condensate 



1 



n = 



1 



(30) 



From a physical viewpoint, the single vortex in a uniform condensate precesses under the influence of its opposite 
image vortex, which explains why the precession rate increases as the vortex moves toward the outer boundary. 

An expansion of A£' (r\) from Eq. 12: ill for small r\ near the center of the condensate shows that the curvature 
changes sign at Q m — 1. This value is the same as the precession frequency for a nearly central vortex in a nonrotating 
condensate; Q m corresponds to the onset of metastability for a central vortex [24| • 

The situation is quite different for the TF condensate, where Eq. (12(:i[l has ln(i?/£) in a term that depends on the 
coordinates of the vortices. Use of Eq. l(2^|l readily yields the corresponding TF precession frequency for a single 
vortex at r% in a nonrotating condensate 
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Mf 



(31) 
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Equation (|31|l differs from Eq. I|30|) because of the "large" factor ln(i?/£), which arises from the radial derivative of 
the energy in the numerator of Eq. (|29|) . Although both 13UI) and (|31|l have similar denominators, the present factor 

(l —rf) 1 reflects the nonuniform TF density p5l l2ri I2?l.l28| (instead of the image, since the TF condensate has no 
image vortices) . Physically, the precession here arises from the combination of the Magnus force and the gradient in 
the density (or equivalently the trap potential in this TF limit) p9ll30ll3i|. 

An expansion of A£' (ri ) from Eq. H26|) for small r\ near the center of the condensate shows that the curvature now 
changes sign at f2 m = ln(i?/£) — i, which is again the same as the precession frequency for a nearly central vortex 
in a nonrotating condensate. This value also corresponds to the onset of metastability, now enhanced by the large 
logarithmic factor. 



E. State of lowest free energy 



In the presence of any small dissipative process, the system seeks the state of lowest free energy for a given external 
rotation by adjusting the radius of the ring of vortices to be r>, determined by the condition 



o= dA£ ' 



dr 3 



(32) 

{rj}=r r 



In this case, Eq. I|29|) shows that = 0, which means that the position of each vortex does not change. Hence all 
vortices are stationary in the co-rotating reference frame. In the laboratory (inertial) frame, the vortices all rotate at 
the same rate as the external rotation fl. Physically this result can be understood when we imagine that the vortex 
core contains some normal component, and the relative motion of the normal component with respect to the frame 
of the rotating trap gives rise to dissipation. 

This equilibrium radius r r is a function of the external rotation f2. For a given vortex configuration, we can evaluate 
the lowest free energy A£^ in (r r , Q) by solving Eq. Comparing these lowest free energies among the different 

vortex states, we can determine the state with lowest energy and the associated angular velocity $l^ B , at which the 
transition from vortex state A to vortex state B should occur. Here we will consider two types of vortex configurations 
for each (uniform/TF) condensate; a ring of vortices without a vortex at the center ('A r ' configuration) and a ring of 
vortices with a vortex at the center ('A-r+Q' configuration). 

For both the uniform and the TF condensates, the vortex-free state is energetically favored for any angular velocity 
below f^ 1 because it has the lowest free energy. As increases, however, the free energy of the one-vortex state 
decreases at a faster rate than that of the vortex-free state. At f^ 1 , the free energy of the two states is the same, and 
the condensate can support a central vortex throughout the interval SI® 1 < ^ < ^l 2 - At higher SI between Oj 2 and 
SI 23 , the state with two vortices has the lowest energy, and so on. 

The two condensates display one interesting difference, however, arising from the different density profiles. In the 
TF approximation, the thermodynamic critical angular velocity SI® 1 (the usual lower critical angular velocity Sl c i) for 
the appearance of the first central vortex is given very generally by 

q01 = A£i = jd 2 r\v 2 n{r) = l jd 2 rr- 2 n(r) 
L\ Jd 2 rn(r) 2 J d 2 rn{r) ' 

where A£\ is the energy per particle for one central vortex and L\ is the angular momentum per particle for the same 
central vortex. Evidently, f^ 1 is proportional to the average value of r~ 2 with the number density as the non-negative 
weight. The TF condensate Utf{t) oc 1 — r 2 has the maximum density at the center; thus the corresponding SI® 1 is 
larger than that for the uniform condensate (the difference leads to a factor ss 2 in the present two-dimensional case) . 

The state with five vortices in a ring plus one at the center (5+® configuration) always has a lower free energy 
than the state with all six vortices in a ring (6 configuration), so the '5' to '6' transition is not favored; instead, '5' 
goes to '5+0'. In fact, this specific '5' to '5+Q' transition was observed in an early experiment on rotating superfluid 
4 He 0. In a uniform unbounded condensate, the vortices form a triangular (Abrikosov) lattice with six-fold rotational 
symmetry. Here the transition from a single ring with an empty center ('A r ') to a ring with a central vortex (W r +®') 
happens when the condensate has five vortices in a ring (N r — 5) instead of six, which reflects the finite size of our 
geometry. 

It is interesting to consider the angular momentum per particle L z , which follows directly from the thermodynamic 
relation 

\ / cq 
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5.1 y 



(a) uniform condensate 



(b) TF condensate 



FIG. 1: Angular momentum per vortex as a function of the external rotation. Here, we assume ln(i?/£) = 8. 



As a function of f2, the envelope of the minimum free energy curves is approximately proportional to — fl 2 , apart from 
small irregularities associated with the small number of vortices, so that L z increases approximately linearly with fi. 
For a given vortex state, the angular momentum increases slowly with the external rotation as the radius of the ring 
of vortices shrinks to satisfy Eq. Since the number of vortices also increases nearly linearly with Q, it is natural 

to consider the ratio L z /N v , which is the angular momentum (per particle) per vortex. This ratio L z /N v is precisely 
1 for a single central vortex. For many vortices, in contrast, the velocity approaches £1 x r (solid-body rotation). In 
our dimensionless units, L z /N v is then the average of r 2 weighted with the particle density n(r), which is | for a 
uniform condensate and ^ for the TF condensate trapped in a harmonic potential. Thus the ratio L z /N v should start 
at 1 for both condensates and decrease slowly toward the appropriate limit with increasing f2. This reduction of the 
angular momentum (per particle) per vortex from 1 was observed in |32j | . 

Figure ^ shows our calculated angular momentum (per particle) per vortex (L z /N v ) as a function of the external 
rotation f2. As expected from the above estimate, the angular momentum per vortex of the TF condensate falls more 
rapidly (typically by ~ 1.4 — 1.6 ~ f) as the number of vortices increases. 



IV. TKACHENKO-LIKE OSCILLATION MODES AND THE ASSOCIATED STABILITY 

Once we have an equilibrium ground state, we can consider small excitations of the ground state by considering 
small variations about the equilibrium. The in-plane, nearly transverse excitation modes of the infinite vortex lattice 
in the uniform Bose condensate are the well-known Tkachenko modes 0. Here we consider the finite-size analog of 
the Tkachenko modes both in the bounded uniform condensate and in the trapped TF condensate. 

Let ctj and (3j be the small variations of the radial distance and the angle, respectively, for the jth vortex in a ring, 



Tj = r r + aj, <j)j = 4>j + (3j, (35) 

where (r r ,<fi®) is the equilibrium position of the jth vortex. In equilibrium N r vortices are symmetrically placed in 
a ring: N r — N v for a state without a vortex at the center and N r = N v — 1 for a state with an extra vortex at 
the center; the state has a A^-fold rotational symmetry. Thus a Fourier decomposition for a and /3 will solve the 
problem: 

(^Sd:)' """' (36) 

where the Fourier index s runs from to N r — 1. 

The vortex at the center requires a special treatment. Since polar coordinates become singular at the origin, it 
is more convenient to use cartesian coordinates (x,y) for the vortex at the center, and polar coordinates (rj,4>j) for 
the rest of the vortices in the ring. In this mixed-coordinate representation with x, y, Tj and <f>j, the Euler-Lagrange 
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equations of motion for the small variations become 

= T^GEid, (37) 

where d = (x, y, a, (3) is a vector. The matrix T2 = diag (2, 2, t, ■ ■ ■ ,t,t, - ■ ■ , t), with t (it is independent of q) defined 
as 



d 2 T 

dr q d<j) q 



(38) 



arises from the time-dependent part T of the Lagrangian. The matrix E2 is the second-order expansion of the free 
energy part AS' 

( 8 2 AS' \ 

<*>«-(svEi),- (3<J) 

and here D = (x, y, r, </>). The antisymmetric matrix G shows how these small- variation parameters couple to form 
the dynamical structure of the system: 






1 




























I Nr 








-lN r 






G = I N \> ^ 
V -I Nr / 

where 7jv r is the N r x N r identity matrix. To apply the above equation (|37|l to the condensate with no central vortex, 
ignore the first two components for x and y in d and set x = y = for the rest of the equations. 

Although these equations look complicated, the Fourier decomposition (|36|l reduces them to a block-diagonal form, 
because A£' includes only two-body terms ~ J* Vj ■ Vj [see Eqs. Q14fl. (I23[l and (|26[l ]. When the Fourier components 
are projected out, the equations for a s and (3 S simplify considerably. For N r = N v — 1 vortices in a ring and one at 
the center, Eq. (|36(l and the orthogonality of the Fourier components decouple the original 27^ x 2N V problem into 
one 6x6 matrix problem, 

^P 6 = M 6 P 6 , (41) 

and N r — 2 sets of 2 x 2 matrix problems (s = 0, 2, . . . , N r — 2), 

^P| = M|P|, (42) 

where Pq = (x,y,ai,aN r -i,Pi,/3N r -i) and P| = (d s ,/5 s ) are 6- and 2-component vectors, respectively. If there is no 
vortex at the center, the original problem corresponds to N r = N v sets of the above 2x2 matrix problems [Eq. H42|) 
with s = 0, 1, . . . , N r — 1]. The matrix elements of M§ and M| arc determined from Eq. 13 7|) (see Appendix iBl) . 

By solving these matrix problems, we can study the stability of the vortex system. For a given equilibrium state, 
we determine the excitation spectrum as a function of the external rotation frequency tt. If a particular excitation 
frequency is real, the corresponding small amplitudes (x,y,a,/3) oscillate around zero, and the state is stable. If any 
excitation frequency is complex with a positive imaginary part, the corresponding mode will grow with time, and the 
state is dynamically unstable. 

As an example, Fig. shows the positive imaginary parts of the excitation spectra for the system of six and seven 
vortices in a ring without a central vortex in the TF condensate. Since the complex excitation frequencies always 
appear in conjugate pairs, the figure shows only the positive imaginary parts. Above a critical angular velocity 
(~ 27.7), a ring of six vortices becomes stable. In contrast, the ring of seven vortices is always unstable because a 
mode with positive imaginary part persists for any Q. 

The Fourier index s is analogous to the wave number in a linear array. For a given N r , the normal modes become 
increasingly unstable as s increases, so that the essential instability arises for short wave lengths. In the case of six 
vortices, s — 3 (the mode with the shortest wave length) is the last mode to become stable. In the seven-vortex state, 
the shortest wave length modes are s = (3 or 4), and these modes remain unstable for all angular velocity, indicating 
a "buckling" transition to a '6+©' configuration, arising from the preference for a more uniform configuration. 
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(a) six vortices in a ring (b) seven vortices in a ring 

FIG. 2: Imaginary parts of the excitation spectrum for a TF condensate showing the normal-mode index s for each curve. 
Here, we assume In (-R/C) = 8. 

A. Uniform condensate 

The dynamical instability of a finite number of vortices in a single ring without a central vortex was studied by 
Havelock 0. There, by solving the set of 2 x 2 matrix problems l|42(l . he concluded that a single ring is unstable when 
the number of vortices exceeds six. Here, we sharpen this result. Although six vortices in a ring are dynamically 
stable, we find that the configuration with '5+©' has lower free energy compared to '6'. 

Figure ES (a) summarizes the conclusions of the free energy calculation and the stability analysis in a uniform 
condensate in a rigid cylinder. Here, the backslashed blocks show the regions where a given configuration has the 
lowest free energy. In the empty region to the left, there exist energetically more stable states with fewer vortices. 
For example, the empty region for the one- vortex state means that the vortex-free state is preferred to the one- vortex 
state. The slashed regions at the left side of the empty regions indicate the absence of a local free energy minimum. 
Hence, a state is metastable, but not globally stable, in the empty region (to the right of the slashed region). For each 
configuration, the narrow solid bar denotes the region where the particular array has at least one complex excitation 
frequency with a positive imaginary part (and is thus dynamically unstable). 

Unlike the ring with no central vortex, the ring with a central vortex is stable even for a ring of eight vortices. 
The dynamical instability appears only for nine vortices in the ring with a central vortex, indicating the existence of 
other more stable states than this '9+®' configuration. Campbell and Ziff 5] studied the stationary vortex patterns 
in a uniform condensate, evaluating the corresponding free energy for numbers of vortices N v — 1 to 30 (they did 
not investigate the small-amplitude normal modes). According to their findings, the stable state with the lowest free 
energy for a system of ten vortices, has '4+4+2' vortices (three concentric rings; beginning from the outer-most ring); 
the '9+®' configuration has slightly higher free energy and is nearly stable, as opposed to stable. 

B. TF condensate 

The TF condensate shows qualitatively the same behavior as the uniform condensate case, with small differences 
[see Fig.Efla)]. 

1. For a single vortex, the dynamics is dominated either by the image vortex (for the uniform condensate) or by 
the nonuniform trap potential (for the TF condensate). For two or more vortices, however, the vortex- vortex 
interactions become increasingly dominant, and the two systems are expected to become more and more similar. 

2. In a transition from '5' to '5+®' for the TF condensate, the region of complex excitation modes (indicated by 
the solid bar) remains almost unchanged, although the number of vortices increases. 

3. Another notable difference is that, for the '8+©' TF configuration, the region of complex eigenmodes (dynam- 
ically unstable) expands significantly, leaving a narrower region of external rotation frequency over which this 
configuration is metastable. 
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(a) uniform condensate 



(b) TF condensate 



FIG. 3: Stability of single-ring states: backslashed blocks denote the region with lowest free energy; empty blocks denote 
metastable states that are local free energy minima but energetically disfavored because other more stable states exist; slashed 
blocks indicate regions where no local free energy minimum exists; the narrow solid bar denotes the region with complex 
excitation frequencies (and thus dynamically unstable). Here, we assume that ln(i?/£) = 8. 



V. DISCUSSION 



The time-dependent variational analysis based on the Lagrangian formalism provides a direct and flexible method 
to study the dynamics of vortices in a Bose-Einstcin condensate. By analyzing the time-independent free energy 
£' , we determine the various ground states and the associated angular velocities for the transitions from one state 
to the next. A more complete dynamical picture emerges when we consider the full Lagrangian C = T — £' ', for it 
provides the explicit time dependence: for example, the vortex precession frequency, the excitation spectrum and the 
dynamical stability. 

For a given configuration of vortices, the critical Q, for the onset of stability is larger for the TF condensate than 
for the uniform condensate, by approximately a factor of two. As a result, a single ring of N r vortices (without and 
with a central vortex) has a smaller radius in the TF condensate than in the uniform condensate, typically, by ~ 0.9. 

The Lagrangian formalism has the notable feature that it apparently incorporates valuable physical information 
that would otherwise emerge only with significant additional effort. For example, the Lagrangian analysis of the 
low-lying normal modes of a condensate (22j used a Gaussian trial function, but the resulting dynamics reproduced 
the known frequencies even in the TF limit where the density profile has a parabolic form. 

Similarly here, the Lagrangian formalism yields the correct precession of a single vortex, even though the assumed 
stream function is equivalent to purely incompressible flow. To understand this question in more detail, note that 
the velocity field assumed in Eq. (|19f) can equivalently be written as v = — V x (z%), which implies that V • v = 0. 
Since we are merely constructing a trial function, this choice is acceptable, even though the correct equation of 
continuity requires V • (nv) = 0. Nevertheless, the Lagrangian approach incorporates sufficient physical information 
(presumably through the energy £') that the resulting precession in Eq. I|31|l agrees with other more lengthy analyses. 

A direct approach to the precession relies on the full equation V • (nv) = v ■ Vn + nV ■ v = 0. The extra term 
involving Vn (namely the spatial variation in the condensate density) produces additional logarithmic contributions to 
the phase [2^, l30l l3l| , which drive the precessional motion. Thus the precession of a single vortex in a TF condensate 
can be considered to arise from the nonuniform trap potential and the resulting nonuniform TF density. It would be 
interesting to improve the trial function by including these additional logarithmic contributions explicitly. Since the 
variational Lagrangian approach is "exact" for the exact solution of the GP equation, such an improved wave function 
should provide an improved description, but it is not clear whether the dynamical precession rate would be affected. 

As studied in this work, the single-ring configuration is applicable only for a small number of vortices (N v < 9 
in practice) in a TF condensate that rotates relatively slowly. The single ring does have one other very interesting 
application, however, when the condensate has a central hole and thus becomes annular. Such behavior is predicted 
to occur in a rapidly rotating trap with an additional quartic confinement, which allows the angular velocity Q to 
exceed the harmonic trap frequency luj_. When the width of the annulus becomes sufficiently narrow, the condensate 
is predicted to have a one-dimensional array of vortices arranged in a single ring |33l l34j . The techniques developed 
here should apply directly to the dynamics and stability of this interesting configuration. 
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APPENDIX A: EVALUATING THE TF FREE ENERGY DIFFERENCE 

With no loss of generality, consider the case with two vortices at ri = (ri,<pi) and r 2 = (r2,<fe)- The stream 
function that describes the velocity field in the condensate is 



X = Xi +X2 = hi|r - n| +ln|r - r 2 | 



(Al) 



In the free energy difference AS' — AS — VLL Z , the angular momentum is straightforward to evaluate and only AS is 
of concern. 



AS = - J d 2 r n(r)v 2 

= X -\ ' d 2 v n(r) (|V Xl | 2 + | V X2 | 2 + 2V Xl ■ V X2 



(A2) 



The first two terms are the contributions of each vortex, and the last term is the interaction between the two 
vortices. For the first term, an integration by parts gives 



- / d 2 rn |Vxi| 2 = \ f d 2 r V • (n Xi V X i) - \ f d 2 r Vn ■ ( Xl V Xl ) - \ \ d 2 r n X iV 2 Xi- 



(A3) 



Excluding the healing length £ neighborhood around the vortex core, the last term in the integral vanishes because 
ri lies outside the domain of the integral [V 2 X i = 27r<5^ 2 - ) (r — n)]. Correspondingly, the surface integral (the first 
term) has two boundaries, a TF radius and a £ neighborhood. With the expansion of the stream function in angular 
harmonics l|21|l . the integral becomes 



\ { d 2 vn\V Xl \ 2 = l -(l-rl) 



- l|) ln(l- ri) 



+ ~ 2 W - 1) 



(A4) 



A similar approach can be used to evaluate the interaction term. 



APPENDIX B: MATRIX ELEMENTS OF THE EQUATIONS OF MOTION 



The orthogonality of the Fourier components simplifies the matrix elements of the equations of motion. 



/ ±AS' iyx \AS\ yy -i\N r E i\N r E 



Ma 



— -AS 

it~ x F 
-it^F 
t- x E 
V t^E 



--AS' 

t- x F 
t- x F 
-it^E 
it- x E 



\N r E 
Ci 


-A 1 




■\N r E 




-A Nr -i 



±N r F 



±N r F 



~i\N r F i\N r F 
Di 
D Nr _! 
-Bi 
-B Nr _xJ 



(Bl) 



and 



Mo 



C s D s 
-A, -B, 



(B2) 
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where 

t 



4 « = 7E A£ / ( r/* rW ' 1 (B3) 
1 w, 

3=1 
1 W " 

C s = -^A^ 5r ..e^-« s , (B5) 

3=1 
1 «r 

F s = -]TA£'^e^-^ s , (B6) 



t 

3=1 



and A£' pq means 



\dpdq J 



cq. 



For a uniform condensate, E = 1 + r 2 , F = — 1 + r 2 , i = 2r; for a TF condensate, E = — (l — 2r 2 + lnr 2 ), 
F = — r (l — 2r~ 2 - lnr 2 ), t = 2r (l - r 2 ). 
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